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Abstract. We study solutions to nonlinear hyperbolic systems with fully nonlinear relaxation 
terms in the limit of, both, infinitely stiff relaxation and arbitrary late time. In this limit, the 
dynamics is governed by effective systems of parabolic type, with possibly degenerate and/or fully 
nonlinear diffusion terms. For this class of problems, we develop here an implicit-explicit method 
based on Runge-Kutta discretization in time, and we use this method in order to investigate several 
examples of interest in compressible fluid dynamics. Importantly, we impose here a realistic stability 
condition on the time-step and we demonstrate that solutions in the hyperbolic-to-parabolic regime 
can be computed numerically with high robustness and accuracy, even in presence of fully nonlinear 
relaxation terms. 
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1. Introduction. We consider nonlinear hyperbolic systems containing fully 
nonlinear relaxation terms when the relaxation is stiff relaxation and for arbitrary 
late time. The dynamics is asymptotically governed by effective systems which are 
parabolic type and may contain possibly degenerate and fully nonlinear diffusion 
terms. For such problems, we introduce implicit-explicit (IMEX) methods which are 
based on Runge-Kutta (R-K) discretization in time, and we apply this method to 
investigate several compressible fluid models. Importantly, we impose here a real- 
istic stability condition on the time-step and we demonstrate that solutions in the 
hyperbolic-to-parabolic regime can be computed numerically with high robustness 
and accuracy, even in presence of fully nonlinear relaxation terms. 

The outline of this paper is as follows. In the rest of this introductory section, 
we review several model problems with linear and nonlinear relaxation, including 
Kawashima-LeFloch's model (TS]. In Section 2, we classify IMEX R-K schemes pre- 
sented in literature and indicate some limitations of these schemes when applied to 
hyperbolic systems with diffusive relaxation. Section 3 is devoted to an introduction 
of partitioned and additive Runge-Kutta schemes for linear relaxation models [SJ [5] . 
In Sections 4 and 5, a particular type of IMEX R-K scheme, called AGSA(3,4,2) (and 
introduced first in [8]), is analyzed and applied to the Euler equations with linear 
friction and to a model coupling the Euler equations with a radiative transfert equa- 
tion (investigated earlier in [5] with a different numerical method). In Section 6, we 
describe and analyze the nonlinear relaxation model from |15j . Concluding remarks 
are finally presented in Section 7. 

1.1. Models of interest. Relaxation effects are an important feature of models 
arising, for instance, in compressible fiuid dynamics and, specifically, we are interested 
here in nonlinear hyperbolic systems with fully nonlinear relaxation terms in the limit 
of, both, infinitely stiff relaxation and arbitrarily late time. In this limit, the dynamics 
of the flow is governed by effective systems of parabolic type, whose diffusion might 
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also be degenerate and/or fully nonlinear. Solving such equations numerically is 
very challenging due to the stiffness of the problem in, both, the convection and the 
relaxation parts. This problem is currently under very active study and the main 
challenge is to design schemes that allow for realistic stability conditions on the time- 
step while ensuring robustness and high-order accuracy. For background as well as 
recent material, we refer to [2l [8l [M l [T3l [T8l [20] . 

The simplest example of interest, which we propose to refer to as the linear 
relaxation model, reads 



ut + v^ = 0, 
vt + h{u)x = -u + q{u) 



2 , , (1-1) 



in which u = u{t,x) and v = v{t,x) are the unknown functions, and b = b{u) and 
q = q{u) are prescribed with b'{u) > 0. This is a nonlinear hyperbolic system with two 
(distinct and real) characteristic speeds, that is, ±y/b'{u)/e. In the stiff relaxation 
limit e — !• 0, the characteristic speeds approach infinity and the behavior of solutions 
to is (formally, at least) governed by the effective system 

ut+q{u)x^b{u)^^, ^^^^ 
V = q{u) - b{u)x, 

which can be derived by a Chapman-Enskog expansion in the parameter e and turns 
out to be a noninear parabolic system — since the diffusion coefficient b'{u) > in 
{b'u)ux)x is positive. Observe that the so-called sub-characteristic condition [23] for 
system (|1.2|) reads e \q'{u)\ < (6'(m))^/^ and puts a restriction on the speed q'{u) of 
the effective equation with respect to the speeds ^y'VJu) / e of the original system: 
clearly, for this model, the sub-characteristic condition is automatically satisfied (for 
e sufficiently small). 

Fully nonlinear relaxation terms also arise (for instance in presence of nonlinear 
friction) and, therefore, in this work we also numerically study the following nonlin- 
ear relaxation model, first introduced by Kawashima and LeFloch [T5] . 

Ut + Vx ~ 0, 

fl 3) 

e^vt + b{u)x = ~\vr-^v + q{u), 

where m > is a real parameter. This, again, is a strictly hyperbolic system of 
balance laws, provided b'{u) > 0. The effective system associated with this model is 
found to be (with a = — 1 + 1/m) 



ut^{\- q{u) + b{u)xr{~q{u) + b{u)x)) ^ 
V = q{u) - b{u)x, 



(1-4) 



which is a fully nonlinear parabolic equation in u. It is natural to distinguish between 
three rather distinct behaviors: 

sub-linear : < m < 1, 

linear : m = 1, (1-5) 
super-linear : m > 1. 

The equation (|1.4p . although parabolic in nature, need not regularize initially discon- 
tinuous initial data, and we may expect jump singularities (in the first-order derivative 
Ux) to remain in the late-time limit. 
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Further models arising in compressible fluid dynamics will also be described later, 
but (jl.ip and (|1.3p provide us two prototypes of particular importance in order to 
develop efficient numerical tools for relaxation problems. 

2. Implicit-Explicit (IMEX) Runge-Kutta schemes. This section is de- 
voted to a brief introduction on IMEX-Runge-Kutta schemes. More detailed descrip- 
tion will be found in several papers on the topic, such as [U l9l [21]. 

Let us consider a Cauchy problem for a system of ordinary differential equations 
of the form 

V' = I{V)+9{V), 2/(0) -yo, ie[0,T] (2.1) 

with y{t) G M^, and /, g two Lipschitz continuous functions — > R^. 

Since Runge-Kutta schemes are one-step methods, it is enough to show how to 
compute the solution after one time step, say from time to time At = h. An 
Implicit-Explicit (MEX) Runge-Kutta scheme has the form 

i — 1 ^1 
Yi=ya + h^dijf{to + cjh, Yj) + h'^aij-g{to + Cjh, K,), 

. 1 
2/1=2/0 + h^bif{to + Cih,Yi) + h^bi-g{to + Cih,Yi). 

i=l 1=1 

where A = (aij), dij = 0, i > i and A = (aij): s x s are lower triangular matrices, 
while c, b,c,b G M." are s-dimensional vectors. IMEX-RK schemes schemes are usually 
represented by a Double Butcher tableau: 



A 



A 



The fact that the implicit scheme is diagonally implicit (DIRK) makes the imple- 
mentation of IMEX simpler, and ensures / is effectively explicitly computed. Order 
conditions for IMEX schemes can be derived by matching Taylor expansion of ex- 
act and numerical solution up to a given order p, just as in the case of standard 
Runge-Kutta schemes. 

In addition to standard order conditions on the two schemes, additional coupling 
conditions may appear, whose number increases very quickly with the order of the 
scheme. However, if c = c and b = b, then there are no additional conditions for 
IMEX-RK up to order p = 3. For more details on the order conditions consult [5]. 

2.1. Classification of IMEX R— K schemes. IMEX R-K schemes present in 
the literature can be classified in three different types characterized by the structure 
of the matrix A = {aij)f of the implicit scheme. Following [3], we will rely on the 
following notions. 

Definition 2.1. An IMEX R-K method is said to be of type A (cf. FH] /) if 
the matrix A G R'^^'^ is invertible. It is said to be of type CK (cf. f^J if the matrix 
A E R"^ can be written in the form 



A 





a A 
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in which the matrix A G invertible. Finally, it is said to be of type 

ARS (cf. m) if it is a specAal case of the type CK with the vector a = 0. 

The schemes CK are very attractive since they allow sonic simplifying assump- 
tions, that make order conditions easier to treat, therefore permitting the construction 
of higher order IMEX R-K schemes. On the other hand, schemes of type A are more 
amenable to a theoretical analysis, since the matrix A of the implicit scheme is in- 
vertible. 

2.2. A Simple Example of first IMEX R K schemes. The simplest IMEX- 
RK are obtained combining explicit and implicit Eulcr method. Two versions are 
possible applied to system (|2.ip . namely 



SP(1,1,1), lU: + (2.2) 



ARS(1,1,1), H]: yi = + h[f[yo) + g[y^)). 
The corresponding Butcher tables are reported below: 



(2.3) 



SP(1,1,1) 









1 


1 




1 


1 



ARS(1,1,1) : 





















1 


1 





1 





1 




1 










1 



2.3. Limitations of IMEX Runge Kutta schemes. In general, Implicit- 
Exphcit (IMEX) Runge-Kutta (R-K) schemes [H 13 H 13 [21] represent a powerful tool 
for the time discretization of stiff systems. However, in the hyperbolic-to-parabolic 
limit under consideration here, the characteristic speeds of the hyperbolic part is 
of order 1/e and standard IMEX Runge-Kutta schemes developed for hyperbolic 
problems with stiff relaxation [5T1 [7| are not very efficient, since the standard CFL 
condition requires the relation Ai = 0{tl^x) between the time mesh-size and the 
space mesh-size. Clearly, in the diffusive regime e << Ax, and we are led to a too 
restrictive condition, since it is expected that, for an explicit method, a parabolic-type 
stability condition At = 0{lS.x^) should suffice. 

Most existing works [13 [13 [20] on asymptotic preserving schemes for nonlinear 
hyperbolic system with diffusive relaxation arc based on a (micro-macro) decompo- 
sition of the hyperbolic Eux into stiff and non-stiff components: the stiff hyperbolic 
component is added to the relaxation term, which is treated imphcitly; in the limit 
of infinite stiffness, such schemes become consistent and exphcit schemes for the dif- 
fusive limit equation [T3 [13 [13 [13 [10] and, therefore, suffer from the usual stability 
restriction At = 0{Ax^). 

Schemes that do avoid the parabolic time-step restriction and provide fully im- 
phcit solvers (in the case of transport equations) have been analyzed in [S], where 
a new formulation of the problem in the case b{u) = u and q{u) = 0, was 

introduced, which was based on the addition of two opposite diffusive terms in the 
stiff limit. One term is added to the hyperbolic component and makes it non-stiff 
(therefore allowing for an explicit treatment), while the other term must be treated 
implicitly. The remaining component of the hyperbolic term is formally treated im- 
plicitly. The resulting scheme is consistent and relaxes to an implicit scheme for 
underlying diffusion limit, therefore avoiding the typical parabolic restriction of pre- 
vious methods. 
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The drawback of the above strategy is that its requires an implicit treatment of 
some components of the hyperbohc flux. Even if such imphcit term can be computed 
by closed formula in most cases, this procedure requires a reformulation of the dis- 
cretization of the whole system. It is thus desirable to construct IMEX R-K schemes 
in which the whole hyperbolic part is treated explicitly. At a first sight this seems an 
impossible task since the characteristic speeds are unbounded. Yet, in [8], the authors 
showed that, provided suitable conditions arc imposed on the IMEX R-K coefficients, 
then one can indeed construct schemes which are 

• fully explicit in the hyperbolic part, and 

• converge to an implicit high-order scheme in the diffusion limit. 

This is this direction that we follow in the present work and further develop in the 
context of, both, the linear relaxation model the fully nonlinear model (|1.3p . and 

several models arising in fiuid dynamics. We exhibit new decompositions that lead us 
to high-order asymptotic preserving methods covering the hyperbolic-to-parabolic 
regime under consideration. A key ingredient is to keep the implicit feature of the 
scheme with respect to certain components that "degenerate" ; for instance, (|1.3|) when 
m > 1 has a degenerate diffusion when ~q{u) + b{u)x is small or even vanishes. 

3. Partitioned and Additive Runge Kutta schemes for the Hnear re- 
laxation model. In developing new IMEX R-K schemes in [8j[6], the system (|l.ip 
was considered and written in the form 



e^vt = -{b{u)x + v - q{u)) 



2.. _ , .J..^^ (3-1) 



A very simple scheme for the numerical treatment of such a problem is based on the 
observation that the first equation does not contain any stiffness, so the first variable 
can be updated in time, and the new value can be replaced in the second equation. 

By using a method of lines approach (MOL), we discretize system (|3.ip in space 
by a uniform mesh {xi}^^-^^ and Ui{t) ?» u{xi,t), Vi{t) w v{xi,t). We obtain a large 
system of ODE's 

Ut=F{U,V), e^Vt=G{U)-V, 

where, in this case, F{U, V) — —VjiV, G{U) = q{U) — 'Djib{U), T>h denoting a discrete 
space derivative operator. In the limit e — >■ the system becomes 

Ut = F{U), where F{U) ^ F{U, G{U)). 

The IMEX Euler scheme (/refAfirst) applied to the system becomes 



yn+l ^ + /\tG{V+^) - Atl/"+\ 



(3.2) 



and the second equation can be explicitly solved for 1^"+^ 

yn+l ^ 



£2 + At 

where we have discretized the interval of integration by a time mesh {inj^i and 

[/" « C/(t"). As e ^ 0, = G([/"+i) and therefore C/"+i = ?7" + AtF([/"), 

i.e. the scheme automatically become an explicit Euler method applied to the space 
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discrctized limit equation. Note that in system (|3.2[) . since the variable U is updated 
in the first equation, it is explicitly computed, although it is formally implicit. 

In the case b(u) = u and q{u) = 0, the method would relax to explicit Euler scheme 
applied to the diffusion equation, thus suffering the usual parabolic CFL stability 
condition At ^ A2:^/2. This approach will be denoted as partitioned approach and 
we will denote the corresponding IMEX R-K schemes, IMEX-I R-K schemes . 

The simplest scheme outlined above is only first order accurate in time. However, 
higher order extensions are possible, by adopting IMEX type time discretization. In 
the previous approach there may be a subtle difficulty when it comes to applications, 
namely it is not clear how to identify the hyperbolic part of the system, i.e. what is 
the term that should be included in the numerical flux. 

For practical applications, it would be very nice to treat the whole term containing 
the flux explicitly, while reserving the implicit treatment only to the source, according 
the following scheme: 



Ut = 














- v/e^ 




[Explicit] 


[Implicit] 



(Additive) 



We call such an approach additive and the corresponding schemes are denoted IMEX- 
E. Such schemes should be easier to apply, since the fluxes retain their original inter- 
pretation. However, this approach seems to be hopeless, due to diverging speeds. 

Let us consider the simple relaxation system with h{u) = u, and q = 0. By 
applying the Implicit-Explicit Euler scheme ARS(1,1,1), (|2.3p we obtain 

=u" + Atv2, 

Now, as e — > 0, one has v"^^ ^'^x- Replacing this expression in the first equation 
gives 

which provides a consistent discretization of the limit equation. 

Performing Fourier stability analysis on this scheme, one obtains a stability con- 
dition that depends on e. As £ one obtains 

At^^ < ^^2 

or 

continuous in space central difference. 

Unfortunately, when trying this approach with other IMEX schemes available in the 
literature, even with the first order SP(1,1,1), one observes a lack of consistency (for 
details, see[8|). 

In this paper, the authors concentrated on developing IMEX R-K schemes of 
type A, since they are easier to analyze with respect to the other types. They started 
the analysis by introducing a property which is important in order to guarantee the 
asymptotic preserving property. IMEX R~K schemes that satisfy this property are 
called globally stiffly accurate, that is the last row of matrices A and A are, respectively, 
equal to iP" and iF . 

Definition 3.1. An IMEX R-K scheme is said to be globally stiffly accurate 



A, b^^e'^A = (0, ...,0, 1)^, c. 
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i.e. the numerical solution coincides with the solution obtained in the last internal 
step of the Runge-Kutta scheme. 

Note that this condition is satisfied by ARS (1,1,1), but not by SP (1,1,1), 

(|2.2p . An IMEX-E R~K scheme applied to this approach will relaxes to an explicit 
scheme for the limit parabolic equation, thus suffering from the classical parabolic 
CFL stability restriction on the time step. 

3.1. Removing parabolic stiffness. The schemes constructed with the ap- 
proachs outlined above will converge to an explicit scheme for the limit diffusion 
equation, and therefore they are subject to the classical parabolic CFL restriction 
At < CAx^. In order to overcome such a restriction we adopt a penalization tech- 
nique, based on adding two opposite terms to the first equation, and treating one 
explicitly and one implicitly. 

Let us consider system (jl.ip and adding and subtracting the same term ^[e)b{u)xx 
on the right hand side we obtain: 

ut= -{v + fj.{£)b{u)x)x + fi{e)biu)xx, 

1 (3.3) 
Vt= — -{b{u)a; + V - q{u)). 

Here /z = ii{e) G [0, 1] is a free parameter such that /Lt(0) — 1. The idea is that, since 
the quantity v + b{u)x is close to q{u) as e — >• 0, the first term on the right hand side 
can be treated explicitly in the first equation, while the term p(u)xX will be treated 
implicitly. This can be done naturally by using an IMEX R-K approach. 

However, if we want the scheme to be accurate also in the cases in which e is not 
too small, then we need to add two main ingredients: 

1. no term should be added when not needed, i.e. for large enough values of 
e, since in such cases the additional terms degrade the accuracy; this can 
be achieved by letting /x(e) decrease as e increases. A possible choice, for 
example, is 

/i(e) = exp(-e^/Ax). 

2. when the stabilizing effect of the dissipation vanishes, i.e. as /.t — > 0, then cen- 
tral differencing is no longer suitable, and one should adopt some upwinding; 
this can be obtained for example by blending central differencing and upwind 
differencing as 

Now from a numerical point of view, we may treat system p.3p in two different 
approaches according to whether the hyperbolic term b{u)x in the second equation 
is treated implicitly or explicitly. The first one, called the BPR approach (For more 
details as well as some rigorous analysis, see [6].) is 



-{v + ^i{€)b{u)r^)x + n{e)b{u)^ 

1 Explicit I 



Implicit 



-{b{u)x + v - q{u)) 



(3.4) 



Implicit 



and the corresponding IMEX R-K schemes are IMEX-I R-K ones in order to remind 
that the term containing b{u)x in the second equation is implicit, in the sense that it 
appears at the new time level. Observe that the term b{u)x + v — q{u) appearing in 
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the second equation is formally treated implicitly, but in practice the term b{u).j. is 
computed at the new time from the first equation, so it can in fact explicitly computed. 
The second approach (cf. [S] for further details) is given by 



-{v + /i(e)6(u)j;)^ 


+ 

_ Explicit 




Explicit 


{v - q{u)) 



XX 



Implicit 



(3.5) 



Implicit 



which is refered to as the BR approach; he corresponding schemes are IMEX-E. 

Observe that when q{u) ^ 0, the system relaxes to a convection-diffusion equa- 
tion. The terms on the left hand side are treated explicitly, while the terms on the 
right-hand side are treated implicitly. In fact as e — > 0, the scheme becomes an IMEX 
R-K scheme for the limit convection-diffusion equation, in which the convection term 
is treated explicitly and the diffusion one is treated implicitly [6l [8] . 

3.2. Additional order conditions. In [HI [5] additional order conditions, called 
algebraic conditions are derived by studying the behavior of IMEX R-K schemes 
applied to the relaxation system p.3|) when e — > 0, in particular when e = 0. These 
algebraic conditions guarantee the correct behavior of the numerical solution in the 
limit e — > and maintain the accuracy in time of the scheme. We obtained such 
algebraic order conditions using the same technique adopted in U, i.e. by comparing 
the Taylor series of the numerical solution with that of the exact solution. 

In [B] these order conditions are derived for IMEX-I R-K schemes and in [5| for 
IMEX-E R-K ones. Due to the number of these order conditions, it is necessary to 
use several internal levels for IMEX R-K schemes in order to obtain a given order. 
Several results and a rigorous analysis about that can be found in [6l [8] . Most nu- 
merical tests are reported in for the BPR approach and in [8] for the BR approach 
and the results are compared with those obtained by other methods available in the 
literature. 

Here we report, as an example, the application of the second order IMEX-E R-K 
scheme developed in [S] and called AGSA(3,4,2). We note that this scheme satisfy all 
the algebraic order conditions derived in |S]. Now we perform several numerical tests, 
considered earlier in [5] with a different numerical strategy. 

4. Euler equation with linear friction. The Eulcr model with friction reads 
ept 4- {pw)x = 0, 

e{pw)t + ipw'^ +p{p))x = -~-pw, 



„2 , „.„^^ _ 1„... (4-1) 



where p > denotes the density and w e R the velocity of a compressible fluid. The 
pressure function p : M+ M+ is assumed sufficiently regular and satisfies p'{p) > 
so that the first-order homogeneus system associated with (|4.1[) is strictly hyperbolic. 
Now using the new variable v = w/e the system becomes 

Pt + ipv)x = 0, 

ipv)t + (pv H = — -pv. 

gZ gZ 

In this case the diffusive regime associated to the Euler model with friction is described 
when e — ^ by the equation 



Pt = [p{p))xx- 



(4.3) 
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Now, from a numerical point of view in order to test our examples we adopt a 
technique similar to the one illustrated in [SJ|H]- Indeed, the more classical techniques 
for instance proposed in [101 1141 [5] provide in the asympotic diffusive limit a scheme 
that converges to an explicit scheme for the parabolic equation. But such schemes 
suffer from the standard CFL restriction At = 0{Ax^) when e — )• 0. Then, in order 
to remove such restriction as suggested in [6l |8] , this technique consists in adding and 
subtracting the same term to the first equation of system ()4.2|) : 



Pt = -{pv + fi{e){p{p)^)^ - fi{e)p{p)^^, 

f ^ / 2^Pip)\ 1 (4-4) 

Here p, is such that /i(e) : M+ [0, 1] and /i(0) = 1. When e is not small there is no 
reason to add and subtract the term {p{p))xx^ therefore /i(e) will be small in such a 
regime, i.e. /i(e) = 0. A possible choice for /i(e) can be 

= (4.5) 
\0, e > Ax. ^ ' 

With this approach, the idea is that as e — 0, the quantity [pv) + p{e)p{p)x ~> 0. 
Therefore such a term can be treated explicitly, while the other term p{p)xx will 
be treated implicitly in the first equation, i.e. we will treat the hyperbolic part of 
the system {{{pv) + fi{e)p{p)x)x, {pv^ + ^^)x)'^ explicitly, while the relaxation one 
{{p{p))xx-,-jiPvY implicitly, respectively. 

We performed our scheme AGSA(3,4,2) introduced in the previous section for the 
numerical experiment considering a parabolic At. i.e. At oc Ax^ by solving system 
(j4.2|) and hyperbohc At, i.e. At oc CFL At by solving (|4.4p . We show only the last 
case since we obtained similar results with the first one. 

We now solve the (|4.2p with the inital data 

'(2,0)^, a; e [1.2,1.8], 
(1,0)-^, otherwise. 



{p,pvr = {)y>^ " (4-6) 



Furthermore we choose the simple pressure law p{p) ~ p^ and e ~ 10^'^. In Figure [SJ 
the numerical solution on a 300 cells grid with Ace = 10~^ obtained by the proposed 
scheme, is compared at time t final = 2.10^e, with a numerical approximation of (j4.3p . 
(or see (4.3) in [5]) computed with 600 cells. 

5. Coupled Euler— radiative transfert model. The second example proposed 
here is a system that degenerates into a system of diffusion equation of dimension 
n > 1 and in order to do this we couple the Euler model (|4.ip with the All model 
proposed in [2] as follows 



ept + {pw)x = 0, 

e{pw)t + {pw'^ +p{p))x = ~-pw + -/, 

eet + {pf)x~0, ' ' (5-1) 

eft + {X ({) e)x = ~-J. 
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Now introducing the new variables v — w/e and / — //e we get the system 
Pt + {pv)x = 0, 

(pv)t + (pv H —)x = — jpw + — /, 

et + [pfh = 0, 



(5.2) 



Here k and a denote positive constants. Furthermore the following pressure law: 
p{p) = Cpp^ with Cp and 77 > 1. The asymptotic diffusive regime (i.e. e — >■ 0) of 
the system (|5.2p is given by 



pip^xx , e-xx 

Pt = + 

K 6k, 

e-xx 



(5.3) 



Now similarly as made for the previous example here in order to remove time step 
parabolic restriction we adding and subtract the same term to the first and third 
equation 

Pt = - (^{pv) + ^ (p{p)x + ^"^^^^ + ^ (pip)xx + ^exx , , 

{pv)t + {pv +—)x = ~^pv+^f, 

et-Hipf) + ^).+ "^, 
1 / \ 



Similarly for this system (|5.4p we perform our scheme considering the hyperbolic At, 
i.e. At oc CFL At. We choose initial data given by 

no n f n J^'^' 2: £ [0.45, 0.55], 

P = 0.2, = 0, / = 0, e = . (5.5) 

I 1, otherwise. 

The parameters of the model are n ~ 2, a ~ 1, e ^ lO^'^, Tfinai = 0.029, Cp = 10^'^ 
and r/ = 2. The numerical solution is computed with 100 cells (Aa; — 10"^) and 
compared to Figure [5] with a reference solution obtained solving ()5.3|) 

In Figures 5.1 and 5.2, we observe that the second-order scheme AGSA(3,4,2) 
developed in [5] captures well the correct behavior of the solutions in the diffusive 
regime where the numerical solutions match the reference solution very well. 

6. Kawashima-LeFloch's nonlinear relsLxation model. 

6.1. Objective. We are now in a position to tackle the nonlinear relaxation 
model (|1.3p which is of central interest in the present paper. Without genuine loss 
of generality in the design of the numerical method, we assume that b{u) = u and 
q{u) ~ and, therefore, we study 

Uf + Vx ~ 0, 

* 1 (6.1) 

e^vt + ux^-\vr-^v, 
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"Reference Solution" 




Fig. 5.1. Comparison between reference solution (full line) obtained by the numerical solution 
or \4--3^ of (4-3) in computed with 600 cells Ax = and A.t oc Ax^ = 5.10^^ versus proposed 
scheme solution (dashed line) ASA(3,4,2) computed with 300 cells and A£ = O.lAa; = 1.10~'^. 
The solutions are compared at time Tj = 2.10^e with e = 10"'^. 



Density 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Radiative Energy 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Fig. 5.2. Comparison between reference solution (full line) obtained by i5.3\ } or the numerical 
solution of (4-6) in ^ computed with 600 cells, Ax = 10"'^ and At oc Ax^ = 5.10~^ versus 
proposed scheme solution ASA(3,4,2) computed with 100 cells and At = O.lAx = I.IO""^. The 
solutions are compared at time Tj = 0.029 with e = 10~^. 

together with the associated effective equation (e — s- 0) 

ut = {\u,\'^u,)^, \v\^-^v^-u,, (6.2) 

with a = — 1 + 1/m. 

The order of convergence of the IMEX scheme apphed to the semihnear system 
(|6.ip for smaU e wih be determined by comparing the numerical solutions to a (the- 
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orctical) limiting solution obtained by solving the effective equations (|6.2p on a very 
fine mesh. 

We are going to treat the nonhncar relaxation term implicitly, which requires to 
solve a nonlinear algebraic equation. An efficient approach is obtained by solving it 
iteratively with the Newton method. Specifically, the nonlinear equation to be solved 
at each s-stage has the form 

F{V) =e^V + CiV"' ~e^C2= 0, 

where Ci and C2 constants. As e ^ 0, we can neglect the term Ve'^ and, therefore, 
we propose the initial approximation 



as a good choice in order to reduce the number of iterations. This is a relevant choice 
in both cases m < 1 and m > 1. 



6.2. Computation of the limit solution. Equation (|6.2p is a nonlinear dif- 
fusion equation. Let us introduce a very efficient method for the numerical solution 
of such an equation. The scheme we propose is stable, linearly implicit, and can 
be designed up to any order of accuracy. We propose to use a technique similar to 
the additive semi- implicit Runge-Kutta methods of Zhong [21]. The idea is to write 
equation as a system 

u'^F{u,u), (6.3) 

where the first argument is treated explicitly, while the second one is treated implicitly. 
Then the semi-implicit R-K method is implemented as follows. First compute the 
stage fluxes for i = 1, . . . , s 

j-i _ i-i 

[/* = w" + At ^ a„ Kj , U, = u" + At ^ a,j , (6.4) 

solve the equation 

K,=F{U*,U, + MauKi), 
and, finally update the numerical solution 



^tY^hK,. (6.5) 



In our case F{u* , u) is given by 

F(u*,u) = (|u*r".)x- 

It can be shown that such semi-implicit R-K scheme is indeed a particular case of 
IMEX R-K scheme [12 . 

Observe that, for the spatial discretization of the limiting equation (|6.2p . we can 
use a centered scheme, such as 
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where \SU\'^_^_ 



U, + iit)-U,{t) 
Ax 



Observe that for the final numerical solution, 

we require, in particular, that u* ^ = it""*"^, which is guaranteed by imposing the 
condition bi = hi for all i. 

Here, as an example, we propose to choose a classical second-order IMEX R-K 
scheme that satisfies the condition 6; = bi for all i, e.g. the implicit-explicit midpoint 
method of type ARS introduced in [1] (to which we refer for further details) . Applying 
such a method, we thus obtain the following (rather simple) scheme which requires 
only two function evaluations 

Un+1 = 2C/2 - Un- 



Observe that U* =Ui= Ur. 

Similar examples is possible to propose for other types of IMEX R-K schemes 
presented in the literature. Then, we consider the numerical solution (|6.6p as a ref- 
erence solution for the limiting equation, and we compare it with the numerical one 
computed from the system (|6.1|) . 



6.3. Comparison with the limit solution. Here we apply first order scheme 
to system (|6.ip with simple second order central differencing for the discretization 
of the space derivative. We choose the first order implicit-explicit Euler scheme of 
type ARS (|2.3p . For our numerical test we integrate the system up to the final time 
T = 1 and use periodic boundary conditions in a; G [— tt, tt] with initial conditions 
u{x, 0) = cos X, w(a;, 0) = sin x and e = 10~^. 

Numerical convergence rate is calculated by the formula 

p = \og2{E^t^|E^t,), 

where i?Ati and i?At2 are the global errors computed with step Ati = 0{/S.x\) where 
Aa;i = 27r/7V and At2 = 0{/\xl) with Ax2 = 7r/iV. 

The nonlinear parabolic equation (j6.2p has regular solutions ifm5^1,i.e. a^O, 
while it develops singularities in the derivatives if0<?7i<l,i.c. a > Q. Two value 
of m are considered here, namely m=l/2(a = l) and m = 2(a = — 1). The profile 
of the solution computed with = 96 points is reported in Fig. 16.31 The time step 
has been chosen to satisfy At — CAx^, with C = 0.025 (right panel with m ~2) and 
C = 1 (left panel with m = 1/2). The scheme seems to converge in both cases. 

Using A^ = 12 • 2P with p = 0,l,...,5, one obtains the convergence table l6Tl The 

convergence is second order for m = 2, and less for to = 0.5, as expected from the 
lack of regularity of the solution. Convergence rate in the latter case improves when 
measured in slightly weaker norms, i.e. and norm. Observe that second-order 
is observed even if the scheme is first-order in time, which is due to the parabolic 
scaling of the time step. 

Let us take a closer look at the case to = 2. By integrating for a longer time, 
some instabilities appear (see Fig. 16. 3|) . The reason of such instabilities is that the 
following equation (|6.2|) may be written in the form 

= ((a + l)|uxr)wxx- (6.7) 
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where the low order term {a + l)\u,j.\" plays the role of the diffusive coefficient. 

It can be shown [5] that the first order Implicit-Explicit Euler scheme applied to 
the linear system (|6.ip suffers from the following stability restriction, independently 
of e 



At 



s; 1, 



which suggests the following condition in the nonlinear case 

At 



(a + l)|M,|" 



Aa 



1. 



This condition is used to determine the optimal time step for m ^ 1, while it shows 
that no time step can guarantee stability near local extrema if m > 1. For this reason, 
it is strictly necessary to use implicit schemes when m > 1. 





(a) m = 0.5 






(b) m = 2 




N 


Error 


Order 


N 


Error 


Order 


12 


1.2942e-01 




12 


7.9684e-01 




24 


5.7400e-02 


1.17 


24 


1.5843e-01 


2.33 


48 


2.2568e-02 


1.34 


48 


3.8728e-02 


2.03 


96 


7.8081e-03 


1.53 


96 


9.3970e-03 


2.04 


192 


2.6057e-03 


1.58 


192 


2.3082e-03 


2.02 


384 


8.2321e-04 


1.66 


384 


5.4599e-04 


2.07 



Table 6.1 

L°° -norms of the relative error and convergence rates of u with C ■ 
C = 0.025 for m = 2. 



1 for case m = 0.5 and 





(a) L^-norm 






(b) L^-norm 




N 


Error 


Order 


N 


Error 


Order 


12 


1.8038e-01 




12 


1.9626e-01 




24 


4.6859e-02 


1.94 


24 


4.0244e-02 


2.28 


48 


1.3354e-02 


1.81 


48 


1.0851e-02 


1.89 


96 


3.575e-03 


1.90 


96 


2.7977e-03 


1.95 


192 


9.3110e-03 


1.94 


192 


7.0519e-04 


1.98 


384 


2.3020e-04 


2.01 


384 


1.6909e-04 


2.06 



Table 6.2 

and -norms of the relative error and convergence rates of u for m = 0.5 with C = 1. 



6.4. Removing parabolic stiflFness for Kawashima-LeFloch model. As 

we have seen, the previous numerical example outlined above converges to an ex- 
plicit scheme for the limit diffusion equation, and therefore subjected to the classical 
parabolic CFL restriction At ^ CAx^. Furthermore in some cases such restriction 
does not allow the computation of the solution, due to the unboundedness of the 
diffusion coefhcient. 

Now, in order to remove the parabolic CFL restriction we adopt a similar tech- 
nique proposed for the linear relaxation model that consists in adding and subtract- 
ing the same term to the first equation. In the nonlinear model equation (jl.3l) . with 
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Fig. 6.1. o numerical solution with TV = 96 cells and solid line the reference limit solution with 
N = 384 cells at time T = 1. The numerical solutions are compared with e = 10~* and At = CAx^. 
On the left-hand side, case m = 0.5 and C = 1. On the right-hand side, case m = 2 and C = 0.025. 





-0.4 -0.2 0.2 0.4 0.6 



Fig. 6.2. On the left-hand side: numerical solution with N = 96 cells at time T = 1.77 for 
m = 2, e = 10~* and C = 0.025. On the right-hand side: zoom of the region where the instability 
appears. 



q{u) ~ 0, b{u) = u, the system becomes 



ut = -{v + iJ.ie)\ux\°'ux)x + fJ.{e){\u^.\°'u^.)j: 



(6.8) 



Here, /J,(e) is such that /.( : [0, 1] and ^(0) = 1. When e is not small there is no 

reason to add and subtract the term fj.{e)uxx, therefore ^(e) will be small in such a 
regime, i.e. /i(e) = 0. 

We note that in both the BR and BPR approach the term {\ux\°'Uj-)a: is treated 
implicitly. In the case of a fully implicit treatment of such term, both approaches can 
be used. Of course the BR approach requires that the Runge-Kutta IMEX is globally 
stiffly accurate O E] . 

If, on the other hand, we want to use the treatment of the implicit term by the 
technique outlined in section 16.21 then we are forced to use BPR approach, since 
the requirement that b — b cannot reconcile with global stiff accuracy. To be more 
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Fig. 6.3. Numerical solution with N = 96 cells at time T = 1.77 for m = 2, e = 10 * and 
At = CAx, with C = 0.25. 



specific, in our case F{y*,y) is given by y = ( " ) ' y* = ( ^* ) > ^^^^ 



F{y*,y) 



V \ V 



It can be shown that such semi-implicit R-K scheme is indeed a particular case of 
IMEX R-K scheme (Pareschi, Russo, private communication). Furthermore, if c = 
Ae = c = Ae, the classical order conditions on the two schemes 

A, 



for order p = 1,2, 3, will guarantee that the semi-implicit R-K has the same order p 
(cf. [5]). In order to show the validity of this technique for the computation of the 
solution for the system (|6.8p . we consider again the convergence test proposed before. 
Here we use a classical second order IMEX R-K scheme of type A with bi = bi, i.e. 
SSP(3,3,2) [21], for the computation of the numerical solution of the system (|6.8p . 

As a reference solution we propose again the numerical solution of the limiting 
equation (|6.3|) and we compare it with the numerical one computed from the system 
(j6.8|) . The L°°-convcrgcncc is second order for m = 2 (cf. Table 16. 3p and a CFL 
hyperbolic condition is used for the time step. 

The following remarks are in order: 

• In the evaluation of the numerical solution for the case m = 2, i.e. a ~ —1, 
we set a tolerance TOL for computing the derivative {\ua: \ + TOL)" in order 
to avoid that the derivatives goes to infinity. We chose in our numerical test 
TOL = 10-12. 

• Observe that by integrating for a longer time with a hyperbolic CFL for the 
case m = 2 the numerical solution decreases rapidity to zero without any 
instabilities (cf. Fig. 16. 3p . 

7. Concluding remarks. In this paper, we have presented a new asymptotic- 
preserving method in order to construct new IMEX Runge-Kutta schemes which are 
especially adaoted to deal with a class of nonlinear hyperbolic systems containing 
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N 


Error 


Order 


12 


1.6921e-01 




24 


4.2166e-02 


2.004 


48 


1.0328e-02 


2.029 


96 


2.5371e-03 


2.025 


192 


6.0394e-04 


2.070 


384 


1.2064C-04 


2.323 


Table 6.3 



L°° -norms of the relative error and convergence rates of u with At = C'Ax and C = 0.06 for 
case m = 2. 



nonlinear diffusive relaxation. As a by-product, effective semi-implicit schemes for the 
limiting (nonlinear) diffusion equations have been proposed. These schemes are able 
to solve the hyperbolic systems without any restriction on the time step which would 
classically be imposed by the stiff source or by the unboundedness of the characteristic 
speeds. In the limit as the relaxation parameter vanishes, the proposed schemes 
relax to implicit schemes for the limit nonlinear convection-diffusion equation, thus 
overcoming the classical parabolic CFL condition in the time step. Although a time 
discretization up to third-order is available [5], we used here second-order schemes, 
since space discretization is limited to second-order accuracy. The construction of 
third-order accurate schemes in space is not a trivial matter in view of the nonlinearity 
and is by itself an interesting field of investigation. 
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